Details of disorder matter in 2D d-wave superconductors 
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We demonstrate that discrepancies between predicted low-energy quasiparticle properties in disor- 
dered 2D d-wave superconductors occur because of the unanticipated importance of disorder model 
details and normal-state particle-hole symmetry. This conclusion follows from numerically exact 
evaluations of the quasiparticle density-of-states predicted by the Bogoliubov-deGennes (BdG) mean 
field equations for both binary alloy and random site energy disorder models. For the realistic case, 
which is best described by a binary alloy model without particle-hole symmetry, we predict density- 
of-states suppression below an energy scale which appears to be correlated with the corresponding 
single-impurity resonance. 
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In metals, the many different analytic and numerical 
techniques used to study disorder and electron localiza- 
tion have led to a satisfyingly consistent understand- 
ing. The same cannot be said for the theory of disor- 
der and quasiparticle localization in high T c supercon- 
ductors, which arc widely believed to be correctly mod- 
elled as 2D systems with d-wave pairing. As pointed out 
some years ago by Nersesyan, Tsvelik and Wenger |l| , the 
standard perturbative approximation for the self-energy, 
the self-consistent T-matrix approximation (SCTMA), 
breaks down in 2D d-wave superconductors as \E\ — > 0, 
even in the dilute impurity limit. In response, a variety 
of non-perturbative approaches have been applied [Q-^), 
and have yielded apparently contradictory results. The 
purpose of this Letter is to demonstrate, by exact numer- 
ical calculation, that these discrepancies occur for the 
most part because, in contrast with the metallic case, 
details of the disorder model are qualitatively important. 
Moreover, for a physically important class of disorder 
models, the seemingly innocent assumption of particle- 
hole symmetric normal state bands leads to non-generic 
results. 

This work focuses on the quasiparticle density-of-states 
(DOS) p(E), which is strongly affected by even small 
impurity concentrations. In pure materials, (i-wave su- 
perconductivity is characterized by a gapless density of 
states p(E) ~ \E\ for \E\ < A (A is the d-wave gap 
amplitude). The SCTMA §-§], predicts a finite DOS 
for disordered materials at E = 0. Non-perturbative ap- 
proaches beyond the SCTMA have variously predicted 
that p(E) vanishes according to universal power 
laws, that p(E) diverges as \E\ — > an d that there 
is a rigorous lower bound on p(E) H. 

At first sight it seems impossible that these results 
could be mutually reconciled. We will undertake to show 
here, however, that most can be understood within a 
single framework, and that they differ primarily because 
of details in the treatment of disorder. It is natural to 
assume that such dramatic discrepancies for what a pri- 
ori appear to be only slightly different physical models 



arise because the (i-wave system is critical. We compare 
models by solving the Bogoliubov-de Gennes equations 
numerically on large finite-size lattices. The model pa- 
rameters we vary include disorder type (see below), dis- 
order strength, normal state particle-hole symmetry, self- 
consistent renormalization of the local order parameter 
by disorder, and the so-called "Dirac cone anisotropy" 
vf/va, where vf is the Fermi velocity and va is the 
velocity of quasiparticles transverse to the nodes. We 
find that binary alloy and random site energy disorder 
models differ qualitatively. For strong scatterers, in par- 
ticular, random site disorder models cannot describe the 
enhancement in the low-energy DOS predicted by Pepin 
and Lee which we reproduce here. In addition, some 
of the present authors have recently shown that self- 
consistent treatment of the order parameter cannot be 
neglected in general. We suggest that an appropri- 
ate model for disorder in the cuprates must involve a 
binary alloy treatment of strongly scattering impurities 
and self-consistency, and predict in this case a power-law 
p(E) ~ E a with disorder-dependent a below an energy 
scale set by the single impurity resonance 0] . 

Method. We consider a mean-field Bogoliubov-de 
Gennes Hamiltonian for electrons hopping on a tight- 
binding square lattice with nearest neighbor hopping ma- 
trix element t, and bond mean field order parameter Ay, 

where U is the impurity potential on site i. Energies 
will be measured in units of the hopping amplitude t 
and lengths in units of the lattice constant. We consider 
both random site energy models, in which the Ui are 
chosen randomly from a distribution P(U), and binary 
alloy models, in which U takes the value Uq on a fraction 
m of the sites and is zero elsewhere. The filling is chosen 
to stabilize a pure d-wave ground state in the absence 
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of disorder, with homogeneous order parameter A& = 
Ao[cos(fc x )-cos(fc a )], where A = |Ei( A »±i _A »±s)- 
In the tight-binding model, vf/va = 2t/A . 

For the d-wave system, Eq. (Q) has been used to calcu- 
late DOS, Qsuperfluid density and T c supres- 
sion numerically It is widely assumed Jl(| that local 
fluctuations in Ay produced by the impurity potential do 
not affect the DOS qualitatively. Recent numerical work 
P,[l2| has suggested otherwise, and we therefore make a 
distinction between self-consistent (SC) solutions of the 



BdG 



where A,- 



equations, >,m^v, 
neighbor pairing interaction Vy, 
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and non-self-consistent 



(NSC) calculations, where Ay has the homogeneous d- 
wave form. With Ay determined, we can calculate the 
DOS from p(E) = L~ 2 Y, a 8(E - E a ), where E a are the 
eigenvalues of 7i, for samples of size L x L. Our numeri- 
cal calculations were performed on systems with L < 45, 
and real periodic and antipcriodic boundary conditions. 

Random site energy vs. binary alloy models. We be- 
gin by discussing the NSC models used in the vast ma- 
jority of earlier work. In field-theoretical approaches it 
is common to assume a Gaussian distribution P(U) = 
(a^/^Tr)^ 1 exp(— U 2 /2<r 2 ) for the disorder potential at 
each site. For technical reasons it is more convenient to 
consider a uniform "box" distribution P{U) — 1/(2W), 
\U\ < W . We have checked that results with the box dis- 
tribution are similar to the Gaussian, with the mapping 
W ~ VScr. Disorder is thus characterized by a single 
parameter, in contrast to the binary alloy model, where 
chemical impurities or vacancies are characterized not 
only by their individual scattering strength Uq, but also 
by their concentration In the normal metal, a corre- 
spondence between W and (ny Uq) can always be found 
such that the random site energy and binary alloy models 
yield similar results. This is no longer true in the super- 
conducting state, because the frequency dependence of 
the superconducting Green's function can lead to midgap 
resonances [ fl3[ , found only in the binary alloy case and 
observed in experiment fl4] , [l5|l . 

Figure |l| shows the DOS of a d-wave superconductor 
at low energies in the presence of box disorder at differ- 
ent W. When W is small, the self-consistent Born limit 
reproduces the exact calculation for the box distribution 
quantitatively. For larger W, the exact calculation shows 
the formation of a "pseudogap" over an energy interval 
\E\ < Ei, where E\ grows rapidly with W. The physics 
of the pseudogap is clearly not captured by the Born limit 
approximation for the box distribution, which predicts a 
finite residual DOS p m [. In Fig. |[ we study the pseudo- 
gap in more detail. With vf = va, arid large disorder, 
we can identify a second, much smaller energy scale E2 
over which p(E) ~ \E\. This regime disappears quickly, 
however, as va is decreased, and we contrast this be- 
haviour with the relatively slow scaling of E\ with va- 
Earlier field theoretical studies || made predictions for 
a linear DOS over an energy scale ~ l/p m f£|. The rapid 
scaling of E 2 in Fig. is consistent with the predicted ex- 




FIG. 1. Density of states for box distributed disorder. 
Main figure: Exact NSC solution of the BdG equations for 
different W. For comparison, the Born approximation for 
W — 5 is also shown. Inset: Born approximation for different 
W. Line types refer to same values of W as in main figure 
(W = not shown). Ao = 2 and p = 1.2. 



ponential dependence of £z, on vf/va- It is also tempting 
to make the connection between the pseudogap edge E\ 
and the much larger predicted scale for weak localiza- 
tion corrections to the DOS 0J§]. In this case, E\ is 
expected to scale as Aop m f for small disorder, and it is 
clear from Fig. |^ that the scaling with Ao holds even for 
large disorder. 

Unitary limit. There is considerable evidence that sim- 
ple defects in the Cu02 planes give rise to local scattering 
centers close to the unitarity limit, indeed that simple de- 
fects in all unconventional superconductors scatter with 
phase shifts close to tt/2, for reasons which are not com- 
pletely understood. For isolated impurities, the signature 
of unitarity is a resonance in the local DOS at E — 0; res- 
onances close to E = have been observed in recent scan- 

H3l for both 
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"native defects" and Zn atoms substituting on the pla- 
nar copper sites in BaSrCaCuO-2212. A point which has 
not been widely appreciated is that the value of the im- 
purity potential Uq which produces a unitary resonance 
is dependent on the band structure; for a perfectly sym- 
metric band the unitary limit corresponds to Uq — > ±00 
|jl3|| , while for an asymmetric band, Uq is a finite value 
dependent on the degree of asymmetry Jl7|-^9[. There 
is a fundamental distinction between the two cases, with 
the first exhibiting perfect particle-hole symmetry on all 
energy scales Q], and the second exhibiting particle-hole 
symmetry on energies \E\ < Aq. Most analytical treat- 
ments of the disorder problem do not distinguish between 
the two, using Aq as a high energy cutoff. We show below 
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FIG. 2. Scaling of the DOS with A for /j = 1.2 and 
= 9. Inset: Logarithmic plots for Ao/2t = 1 (circles) 
and Ao/2t = 0.3 (diamonds; data rescaled for clarity). Lines 
are guides to the eye indicating linearity. 



that in the many-impurity problem erroneous conclusions 
may be drawn as a result. 

In the unitary limit, the SCTMA predicts that a 
plateau forms in the DOS over an energy interval \E\ < 7. 
Recent non-perturbative calculations by Pepin and Lee 
P| found that unitary scatterers produce a divergent 
DOS p{E) ~ l/\E\\vi{\E\/l±o) as \E\ -> 0. This fea- 
ture was not found in recent numerical work by some 
of the current authors (9j, who considered only tight- 
binding bands with \i ^ 0. In Fig. [| we show that a 
divergent DOS occurs only in models with particle-hole 
symmetry at all energy scales, for — and Uq 1 = in 
our case. The effects of breaking particle-hole symmetry 
are illustrated in Fig. |[ For /1 = and Uq 1 7^ 0, the 
divergent peak splits, and moves away from E = as 
IC/q -1 ! grows. This is qualitativel y s imilar to what hap- 
pens in the single impurity limit |l3|j , although the peak 
splitting occurs more rapidly in the bulk disordered case. 
An alternative means of breaking particle-hole symmetry 
is to let /i 7^ 0, in which case the peak structure rapidly 
disappears (Inset, Fig. ||). This differs from the single 
impurity limit where there always exists some finite Uq 
at which a zero energy divergence occurs. A remnant of 
the isolated impurity resonance remains, however, as a 
broad accumulation of states at low energies. The extra 
spectral weight is not reproduced by SCTMA calcula- 
tions, and it gradually vanishes as we move further away 
from perfect symmetry. Evidently the low-energy effec- 
tive action is unstable both to deviations from unitarity 
and to deviations from full band particle-hole symmetry. 

Local order parameter suppression. In the above dis- 
cussion, we have concentrated on the artificial case in 
which the order parameter was not allowed to vary spa- 
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FIG. 3. NSC calculation of DOS near the unitary limit for 
A = 2. Main figure: /j, = 0. Inset: fi = 0.2 with Uq 1 = 0.001 
(solid), 0.01 (dotted), 0.02 (dashed). 



tially, with the purpose of understanding a set of dis- 
parate theoretical results obtained under this assump- 
tion. As discussed in Ref. H, allowing Ay to respond 
sclf-consistently to the impurity potential introduces a 
new source of scattering in the off-diagonal channel, 
which ultimately leads to a suppression of the DOS at 
low energy. This effect highlights the complex nature of 
multi-impurity scattering resonances, and is opposite to 
the expectations of the naive "Swiss cheese" picture, in 
which the impurity simply produces a small region of nor- 
mal metal. In Fig. ^ we show the effect of self-consistency 
in the unitary limit. In the case of a symmetric band, 
self-consistency moves the resonance towards the Fermi 
level (consistent with what is seen in the single impurity 
case [ po|Jl9[| ) and also suppresses the overall low energy 
DOS. This kind of DOS suppression is also seen in SC so- 
lutions with box distributed disorder [Fig. §(b)]. In the 
binary alloy model, however, it is possible to make an 
empirical connection between the energy scale for DOS 
suppression and the energy of the single impurity reso- 
nance. The energy scales are not equal, but clearly scale 
together M . One interesting implication is that for scat- 
terers sufficiently close to the unitary limit, the pseudo- 
gap will be unobservable, as illustrated in Fig. [|(c), but 
that the same impurity doped into a different host may 
cause the pseudogap to open [Fig. f|(d)]. We speculate 
that this may be the case with Zn, which is clearly a uni- 
tary scatterer in BSCCO |l4(] , but appears to produce a 
depression in the DOS at small energies in LSCO, as seen 
in recent specific heat experiments [ pl| . 

Conclusions. We have shown that a number of dif- 
ferent approaches to the d-wave disorder problem which 
produce apparently contradictory results can be under- 
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FIG. 4. Exact solutions of BdG equations with (SC) and 
without (NSC) order parameter self-consistency, (a) Binary 
alloy, symmetric band, near unitarity; Ao = 2, Uo = 100, 
m = 0.05, ju = 0. (b) Box distribution; A = 0.39, W = 2, 
[i = 1.2 (c) Realistic binary alloy model; Ao = 0.39, Uo ~ 5, 
m = 0.04, fx = 0.6. (d) As in (c) but with fx = 1.2. 



stood within a single framework when the appropriate 
symmetries of the Hamiltonian, and in particular of the 
particular realization of disorder, are accounted for. The 
most important question we hope to settle here is which 
of the preceding results, if any, are of relevance to ex- 
periment. We remind the reader that experiments which 
probe the DOS most directly are consistent with the ex- 
istence of a constant DOS at the Fermi energy |p2f . We 
have shown that a true constant DOS cannot be under- 
stood in 2D (i-wave superconductors with any of the dis- 
order models discussed here [^3|, and we suggest several 
possible reasons for the discrepancy. The first possibility 
is that experiments are unable to access the pseudogap 
regime in the optimally doped materials because of na- 
tive near-unitarity defects |Q. In this case, the only ef- 
fects of weak localization on the DOS would be the weak 
nonmonotonicity shown in Fig. [|(c),(d), which is remi- 

|§ The second 
lird dimension 



niscent of effects seen in the specific heat. | 
possibility is that weak coupling to the t 
destroys the DOS anomalies described above; we expect 
on general grounds that the influence of the crossed dia- 
grams identified by Nersesyan et al. Q| will become neg- 
ligible in the dilute limit in 3D, and the validity of the 
SCTMA will be restored. In this context we note that 
almost all the experiments indicating finite residual DOS 
in the cuprates have been performed on YBCO, the most 
3D of the cuprate materials. We hope our work will serve 
as an incentive to examine the low-energy properties of 
disordered 2D materials like BSCCO-2212. Finally, we 



mention the possibility that many-body effects beyond 
the BCS approximation play an important role at low 
energies. There is some speculation that this might be 
the case in the underdoped cuprates, and the formulation 
of a relevant theory is an interesting but difficult problem 
which must be left to future research. 

Note added: In the final stages of preparation of this 
manuscript, we received a preprint from Zhu et al. ]25j] in 
which similar results for the unitarity limit of the sym- 
metric band were obtained. 
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